SEISMIC HALOS AROUND ACTIVE REGIONS: AN MHD THEORY 
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ABSTRACT 

Comprehending the manner in which magnetic fields affect propagating waves is a first step 
toward constructing accurate helioseismic models of active region sub-surface structure and dy- 
namics. Here, we present a numerical method to compute the linear interaction of waves with 
magnetic fields embedded in a solar-like stratified background. The ideal Magneto-Hydrodynamic 
(MHD) equations are solved in a 3-dimensional box that straddles the solar photosphere, extend- 
ing from 35 Mm within to 1.2 Mm into the atmosphere. One of the challenges in performing these 
simulations involves generating a Magneto-Hydro-Static (MHS) state wherein the stratification 
assumes horizontal inhomogeneity in addition to the strong vertical stratification associated with 
the near-surface layers. Keeping in mind that the aim of this effort is to understand and charac- 
terize linear MHD interactions, we discuss a means of computing statically consistent background 
states. Power maps computed from simulations of waves interacting with thick flux tubes of peak 
photospheric field strengths 600 G and 3000 G are presented. Strong modal power reduction in 
the 'umbral' regions of the flux tube enveloped by a halo of increased wave power are seen in 
the simulations with the thick flux tubes. These enhancements are also seen in Doppler velocity 
power maps of active regions observed in the Sun, leading us to propose that the halo has MHD 
underpinnings. 

Subject headings: Sun: helioseismology — Sun: interior — Sun: oscillations — waves — hydrodynamics 



1. INTRODUCTION 

The complexity of the solar background state, subtleties in the dynamics of wave propagation in the 
near-surface layers, and the inherently anisotropic, tensorial nature of magnetic fields disadvantage analyt- 
ically driven MHD studies. There have been many theoretical efforts to model MHD interactions in flux 
concentrations but have proven to be somewhat restrictive in the scope of problems addressed given the ef- 
fort required to construct these models. In this regard, numerica l forward modeling of wave propagation (e.g. 



Hanasoge et al.l2006ll2007al:ICameron Gizon. fe Daiffallan1l2007tlShelvag et al.|[2007tlKhomenko. Collados. fc Felipe 



2007t|Parchevskv fc Kosovichevl2007a ) has been relatively successful at making sense of the sometimes highly 
counter-intuitive wave phenomena observed in the Sun. 

Accurately deconstructing the sub-surface structure and dyna mics of active regions is a difficult task . 
Since the development of methods of time-distance helioseismology ( Duvall et al. 19931 Gizon fc Birch 2005) 



and the subsequent inve stigation s into the n ature of the sunspot underbelly (jDuvall et al . 199G: 



1997 



Kosovichev fe Duvall 



Couvidat. Birch, fc Kosovichevl200tl ). there have been several arguments attempting either to establish 



the s i gnificance of MHD in teractions in sunsp ot structure and dynamics inversions (e.g. Lindsev & Braun 



2005 



Schunker et al. 2003h or to the contrary ( Zhao fc Kosovichev 2006h . Recent theories ( Braun fc Birch 
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20061 ) argue that most of the observed wave phase shifts in sunspot regions occur in a thin sub-photospheric 
region of 1 Mm depth, where magnetic field effects are putatively the largest. The implication is that the 
ca usative mechanisms behi nd observed wave phase shifts may have been misidentified, a conclusion echoed 
by lHanasoge et al.l (|2007bh who demonstrate that wave source suppression due to convect i ve bl ocking in 
sunspots can also participate in the task of creating time shifts (also see, iGizon fc Birchl 120021 ) . More- 
over, wave phase sh ifts inferred in regi ons of strong magn etic fields from Michelson Dopplcr Imager (MDI; 
Scherrer et al.lll995l ) observations (e.g. iDuvall et al.lll996l ) are difficult to interpret because of substantial 
changes in the line formation height due to profound alterations in the thermal structure of the underlying 
plasma. On the positive side, the prevalence of computing resources and numerical methodology now af- 
ford us the ability to conduct investigations that may not have been possible a decade ago. Developing an 
interaction theory of waves and magnetic fields will allow more consistent studies of sunspot structure and 
dynamics. 



The reduction in acoustic oscillation power in sunspot regions (e.g. iLites et al.lll982l) h as been the sub- 
ject of extensive observations with seve r al theories put fo r th to explain this phenomenon (e.g. lHindman et al 



19971 : iParchevskv fc Kosovichevll2007bf ) . iHindman et al.l f| 19971 ) have discussed several plausible mechanisms 
that may be contributing to the power reduction but the participatory extents are as yet unknown. On a 
related issue, a number of st udies have focused on placing observational constraints on the de gree of wave 
absorption in sunspots (e.g. iBraun et al.l Il987t iBogdan et all Il993t iBraunl Il995t ICallvl Il995h . The tech- 
nique discussed here provides an independent manner of investigating all these issues. Acoustic or seismic 



enhancements (or halos as they are termed in t his paper) are ubi c 
sity observations, encircl i ng active regions ( e.g. 



Hindman fc Brown 1998: Doneaetal 



Braun et al 



20001 : iNagashima et al 



992; 



uitously seen in both velocity and inten- 



Brown et al 



20071 ). Some (|Brown et al 



1992 



Ba 



199- 



thasar et al 
19921: toonea et al 



200* 



0) have speculated that they originate from enhanced source activity in the vicinity of the active region 



In this paper, we present power maps from simulations of waves interacting with moderate to strong magnetic 
fields; acoustic halos are clearly seen in these images, implicating an MHD based mechanism. 

On a very different scale but of equal importance are small magnetic elements and thin flux tubes. 
The dynami cal emergence a nd disappearance of these flux tubes provides us insigh ts into the photospheric 
dynam o (e.g. lCattaneolll999i ). In a bid to understand the structure of these flux tubes. lDuvall. Birch, fc Gizon 
(|2006l ) analyzed MDI observations of thousands of independent small magnetic elements, thereby developing 
a highly resolved statistical picture of the associated wave scattering. Understanding the nature of the 
interaction between thin flux tubes and waves may allow us to recover details of the flux tu be structure from 
the scattering information. Forward models of wave interactions with thin flux tu bes (e.g. IBogdan fc Callv 
1995t IBogdan et a.1.1 Il996t IGizon. Hanasoge. fc Birchl 120061 ; lHanasoge et al.l 12007a ) can then be constructed 
in order to place restrictions on the subsurf ace magnetic field di stribution. Models of this sort can be used 
in theoretical studies of flux emergence (e.g. ICheung et al.ll2006l ). 



In this regard, a first step is to devise a sufficiently general manner of computing wave propaga- 
tion in a magnetized plasma. The linearized ideal MHD equations provide a reasonable starting point, 
since MHD o scil lations in the ph o tosph er e and below are governed by pre dominantly l inear physics (e.g. 
Bogdanllioool ). ICallv fc Bogdar] |l997l ). iRosenthal fc Julienl |2000l ). and ICallvl (|2000h performed MHD 
simulations in two dimensions to st udy rates of mode absorption in magnetic flux tubes. Subsequently, 
Cameron. Gizon. fc Daiffallahl (|2007l ) developed and validated numerical techniques to perform 3D linear 
MHD computations with a focus on recovering the magnetic field distribution based on wave scattering 
measurements. The assumption of linear wave propagation and time stationarity of the background state 
are common threads between this work and that of the above-cited authors. 
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High-order numerical accuracy is a minimum requirement for computational work. The linear calculation 
discussed here does not face the same restrictions as would a non-linear counterpart, where the presence of 
shocks makes it quite difficult to raise the order of the numerical scheme without introducing instabilities. 
We discuss the methods employed to spatio-temporally evolve solutions of the ideal MHD equations in §21 
Subsequently, an empirical method to generate stable MHS states is introduced in Sj3j with an illustration 
of one such state: a flux tube with peak photospheric field strength 600 G. Results of wave simulations with 
some flux tubes, specifically the phenomena of wave power reduction and enhancement are discussed in §3.11 
Finally, we summarize and conclude in §21 



2. COMPUTATIONAL METHOD 



(2 



Similar to the forward models of the solar wave field developed in lHanasoge et al.l (|2006f ) and lHanasoge et al 



200IaO, we start by linearizing and modifying the ideal MHD equations in the following manner: 

d t p = -v-( Po v)-r P , (i) 

ftv = - — Vp - -%e 7 + ^-[(VxBo)xB + (VxB)xBo] + S - IV, (2) 

d t p = -v • Vp - p a c 2 V-v - Tp, (3) 

9 t B = CVx (vxB ) — TB (4) 



V • B = (5) 

where p denotes density (unless stated otherwise, the subscript '0' indicates a time-stationary background 
quantity while un-subscripted terms fluctuate), p pressure, B = (B Xl B y , B z ) the magnetic field, v = 
(v x ,Vy,v z ) is the vector velocity, g — g(z) is gravity with direction vector — e 2 , c = c(x,y,z) is the sound 
speed, r = T(x. y,z) > is a damping sponge that enhances wave absorption at all horizontal and vertical 
boundaries (see Figure [l}, £(z) a Lorentz force 'controller' (Robert Cameron, private communication 2007; 
Robert Stein, private communication 2007), and S is the source term. The controller term £ (see Figure^) is 
such that it is constant (=1) over most of the interior but decays rapidly with height above the photosphere. 
Note that £ is also present in equation ^ - as the influence of the magnetic field on the fluid decreases 
(Eq. [5]), so must the effect of the fluid on the magnetic field. For further discussion on the reasoning behind 
this term, see £)2.1I 

We employ a Cartesian coordinate system (x,y,z) with e z denoting the unit vector along the vertical 
or z axis and t, time. Because of the presence of a spatially varying magnetic structure, the background 
pressure, density, and sound speed adopt a full three-dimensional spatial dependence. In sequential order, 
equations |T]) through ([3]) enforce mass, momentum, and energy conservation respectively, while equation ^ 
is the induction equation. Equation |(SJ) assures us that magnetic monopoles do not exist. In interior regions 
of the computational box (away from the boundaries), solutions to the above equations are adiabatic since 
the damping terms decay to zero here. The source term S, is a spatio-temporally varying function, the 
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structu re of which has been discussed in some detail in lHanasoge fc Duvalll (|2007l ) and lHanasoge et al 
(|2007al ). Essentially, it is a phenomenological model for the multiple source wave excitation picture that 
i s observed (inferred p erhaps) in the Sun. The background ver tical stratification is an empirically derived 
( Hanasoge et al.l l2006). convectively stabilized form of model S ( Christensen-Dalsgaard et al.lll99 



The base hydrodynamic method remains unchang ed fro m Hanasoge et al.l (|2007al ); spatial derivatives are 
calculated using sixth-order compact finite differences (|Lelelll992h and time evoluti on is achieved throu gh the 
repeated application of an optimized second-order five-stage Runge-Kutta scheme ( Berland et alj|2006l) . The 
temporal order of accuracy is dropped because the time step (2 seconds) is much smaller than the period of 
the waves studied here. The boundaries are lined with damping sponges in order to absorb (damp) outgoing 
waves (Figure[T]). This is to prevent any scattered waves from re-entering the computational domain as would 
be the case with periodic boundaries. Our attempt to extend the base scheme to compute the magnetic 
field terms in equations (|2|) and (j4|) was successful. All derivatives, including the magnetic field terms, are 
estimated using sixth-order compact finite differences, thus maintaining a high order of spatial accuracy. It 
was observed that the V • B term was of a low magnitude, < 10 -7 per pixel, and therefore harmless (e.g. 
T6t,hll200fi I Abbettl 120071 ) . Moreover, the presence of the damping term TB ensures that V • B is forced to 
decay in the damping sponge layers. Validation in one and two dimensions of the essential numerical method 
(i.e. without the £ or T terms) is discussed in Appendix lAl 




Fig. 1. — The function r(x, y, z) of equations (JTJ - Q. In panel a, the value of F at a location far away from 
the side boundaries is plotted as a function of z. The vertical boundaries of the computational box are at 
z = 0.95, 1.002./?©; although not shown here, the behavior of T at the lower boundary is qualitatively similar 
to that in panel a. In panel b, the variation of T with the horizontal coordinate x at a location far away from 
the vertical and y boundaries, with x — serving as one of the side boundaries. As waves approach within 
20 Mm of the horizontal and/or 1 Mm of the vertical boundaries, they start to experience strong damping 
by the —T term. 
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2.1. Lorentz force controller 



As stated in ^ £ retains the value 1 in the interior and decreases with height above the photosphere 
(Figure^). It attempts to achieve a two- fold purpose: (I) a reduction in the Lorentz force with increasing 
altitude above the photosphere and (II) prevent the onset of negative pressure effects. The mean hydro- 
dynamic pressure and density in Sun drop exponentially with height in the atmosphere that immediately 
overlays the photosphere. In calculations of MHS states (§E]), it was nearly impossible to prevent complete 
pressure and density evacuation in the interiors of flux tubes of large magnitude field strengths (1500 Gauss 
and more - sa dly nowhere close to t he umbral field strengths of up to 6100 G that have been observed 
in sunspots by Livingston et al. 2006h . Moreover, the equilibrium horizontal pressure distribution takes on 
strange forms, with the pressure at the center of flux tube attaining larger values than the ambient, when the 
flux tube radius is forced to increase faster than the corresponding potential field configuration. In the Sun, 
the presence of magnetic field everywhere and the phen omenon of flux tube merging in the atmosphere (e.g. 
Pneuman. Solanki. fe Stenflo 1986 : Bogdan et al. 1996h help reduce large gradients in the magnetic field, 
thereby preventing complete evacuation in active regions and sunspots while not requiring the flux tubes to 
flare out too rapidly. We attempt to simulate this (criterion I) through the £(z) term. However, since the 
equilibrium structure of an active region is as yet unknown, it is not possible to determine how realistic a 
chosen functional form of £ is. 

To determine the impact of £(z) on the wave field, we simulate the interaction of a wave packet with a 
relatively weak flux tube (~ 100 Gauss at the photospheric level) in a solar-like stratified medium. Three 
simulations are performed, a quiet simulation ('q') without any magnetic field and two ('c' and 'd') with 
different functional forms of £, one form of £ decaying more rapidly with height than the other (shown 
in panel a of Figure [2]) ■ The initial condition for all simulations was chosen to be a Gabor wavelet shaped 
disturbance in v z localized at (x, z) = (—30, —0.2) Mm, at all y. At approximately the instant when the wave 
packet reaches the center of the flux tube (located at (x,y) = (0,0) Mm), we display snapshots of v\ — 



and 



in panels c and d respectively, where the superscripts refer to the simulation index (q,c, or d). 



In the presence of a linear scatterer, one may view the velocity field as being associated with an incident 
and scattered wave; in this situation, is the incident wave velocity, while the scattered wave velocities 
are described by the differences v c z ' d — v^. It is clear from panels c and d of Figure [2] that the extent of 
scatter in the simulation where £ decays higher up in the atmosphere (c) is greater by an order of magnitude 
than (d). Perhaps mode conversion, w hich has been theor e tically shown to become significant when th e 
plasma-/? starts to drop, is at play (e.g. iBogdan et al.lll996l : [Callv fc Bogdan] Il997l : iGrouch fc Calrvll2003h . 
It may also be that the magnetic field changes the v z eigenfunction more significantly in one case than the 
other. Essentially, this experiment tells us that capturing wave interactions in an active regions is somewhat 
sensitively dependent on the choice of the £ function, or in other words, on the atmospheric magnetic field 
distribution in the vicinity. It underlines the neccessity of viewing this effort as more qualitative than 
quantitative, since conclusions of the latter sort require exploring a formidable parameter space. 



MHS STATES 



Generating MHS states in stratified media can be a non-trivial task (e.g. iPneuman. Solanki. fc Stenflo 



1986 



Pizzo 



1990; 



Belien et al 



2002 



Khomenko. Collados. fc Felipdl2007l ). Fully consistent approaches that 
involve relaxing the MHD equations to low-energy equilibria arc difficult to implement. Moreover, such 
calculations are beyond the scope of this effort; we are interested less in the MHS state itself than in 
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Fig. 2. — Figure showing the dependence of wave scattering on ((z) (Eq. [2], [4]). Panel a shows functional 
forms of £, used in two different MHD wave test simulations, termed c and d. The initial condition, same 
for both simulations, is a plane wave packet localized at x = —30 Mm, z = —0.2 Mm. Panel b shows 
fluctuations in b z , which arise due to the interaction of the wave with the flux tube. Panels c and d display 
the instantaneous normalized vertical velocity (v z (x,y,z = 0.2 Mm, t = 31 min), units are arbitrary) of the 
scattered waves extracted at a height of 200 km above the photosphere from simulations c and d respectively. 
An order of magnitude difference is seen between the two cases, indicating that the results are somewhat 
sensitively dependent on the chosen form of 
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thc manner in which waves interact with them. We invoke the Schliiter fc Temesvarv ( 19581 ) self-similar 
magnetic field geometry and ignore both radiative transfer effects and the satisfaction of the equation of 
state. We also remind the reader that the background stratification has been altered to prevent the onset of 
unc ontrolled linear growth of conv ective instabilities, thus changing the opacities in a non-physical manner. 
The ISchliiter fc Temesvarvl (|1958T ) approximation tells us that making the following c hoices for the radial 



and v ertical magnetic field B r and B z assures us of the satisfaction of equation (O (e.g 
2005h : 

B z = MiP(z)e- r2 ^ z) , 



Schiissler fc Rempel 



B r = -M-i>'e- 



(G) 
(7) 



with tp' = dtp/dz. The above equations © and |(7J) are in cylindrical geometry; r, z refer to the horizontal 
radial and vertical coordinates with r = coinciding to the center of the flux tube, M a term that controls 
the magnitude of the magnetic field and hence the flux (= ttM), and *p(z), the horizontal extent of the 
flux tube and the rate at which the flux tube spreads with altitude. The zeroth-order MHS equations in 
cylindrical coordinates, obtained upon dropping the time and azimuthal dependencies in equation ((2]), reduce 
to: 

= -d r p + [d z B r - d r B z ] , (8) 



along the horizontal (r) direction and in the vertical (z) direction, 



= -d z p-(^[d z B r 



d r B z 



pg- 



Equation 



is integrated from r = 0, oo to obtain the following equation 

M 2 C 



Pc(z) = Poo (2) + 



4tt 



16 V 



8 W 2 



(9) 



(10) 



where p c {z) is the pressure along the axis (centerline) of the flux tube and Poo(z) is the hydrostatic pressure 
far away from the magnetic region. The horizontal pressure distribution at a given z can now be computed 
by integrating equation (|5|) from the center outward: 



p(r',z) =p c (z) + 



f dvB z [d z B r - d r B z ] ■ 
Jo 



(11) 



thus the entire pressure field can be recovered through this procedure. Simplifying equation ([5]), we can 
obtain the density field from the pressure distribution: 



p{r, z) = — (d zP + [d z B r - d r B z 
g \ 4tt 



(12) 



Therefore, upon specifying parameters M and ^{z) in equations ^ and ijTj). one can obtain a self-consistent 
MHS solution that satisfies the criteria of V • B = and magneto-hydro-static balance. One must be careful 
however to ensure positive pressure and density in equations (fTT|) and (TT2|) at all points in the computational 
domain. In Figure El we show an example of a flux tube which attains a peak strength of 600 G at the 
photospheric level; the inclination of the field at distances away from the center is also shown. 
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Fig. 3. — An example of a flux tube generated according to the recipe of S}3] Panel a shows the field strength, 
\B\ = [B 2 + B 2 ] 1 / 2 ; panel b shows the field inclination &taii(B r /B z ). Perpendicular to the contour lines are 
spokes that point in the direction of the downhill gradient. The pressure and density remain positive over 
the entire domain. Field strength magnitudes in Gauss and inclination angles in degrees are indicated along 
the contour lines. 
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Fig. 4. — Changes in acoustic power caused due to MHD interactions (600G case). Two simulations were 
performed: a quiet ('q') and an MHD calculation ('mag') with the flux tube of Figure [3] embedded in the 
computational domain. Panel a shows the RMS differences of the total velocity, u to t = Jv* + Vy + 
between the quiet and magnetic simulations. Panels b and c display the RMS differences seen in v x and 
v z while d is the RMS of the vertical magnetic field fluctuations, b z . Each power difference is normalized 
by the mean value of the RMS of the corresponding quantity (i.e. v x ,v z , or i>tot) derived from the quiet 
simulation. With increasing radius, the contours in panels b and d show locations where the field inclination 
is [20°, 30°, 40°, 50°]. The field strengths at these contours are [560, 495, 397, 240] G respectively. Outside 
a region of substantial decrease in wave oscillation amplitudes, a halo corresponding to an increase in the 
RMS is seen. 

3.1. Seismic Power Deficits and Halos 

Theoretical expectations dictate a decrease in modal power in magnetic regions due to mode absorption 
and MHD-wave coupling. Using identical realizations of the source function S (Eq. [2]), we perform two 
simulations: a 'quiet' run and an MHD counterpart ('mag') with the flux tube of Figure [3] embedded at 
the center of a computational box of size 100 x 100 x 35 Mm 3 via the computational method of S}2] In 
Figure 21 we show the difference in time-averaged Root Mean Square (RMS) wave power between a quiet 
simulation and its magnetic counterpart, normalized by the mean value of the RMS of the quiet case. Both 
runs were twelve hours long. Because the realizations are identical, the MHD interactions are the dominant 
component of the quantity RMS mag - RMS q . It is interesting to note that depending on the variable of study, 
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the simulations predict strong variations in the nature and degree of change in wave power. For example, 
the RMS differences in the total velocity, v to t - 



of, + + v\ show the presence of a large reduction in 



wave power surrounded by an intense halo, whereas the RMS decrease as seen in v z is systematically weaker 
and an almost invisible halo. The panels b and d show contours of increasing radii corresponding to field 
inclinations of [20°, 30°, 40°, 50°]. The halo is seen at inclinations of 50° and higher, while a strong reduction 
in wave power is observed at smaller angles. Also, the robustness of the halo was ensured by verifying 
its reappearance in a sim ulation using an alterna te numerical method, namely a second-order Constrained- 
Transport technique fCT; lEvans fe Hawlevlll988l ). 

In order to study these effects further, we computed power maps in four different frequency bandpasses, 
2 - 3, 3 - 4, 4 - 5, and 5-6 mHz. Different components of the velocity were used in the calculations, 

u z ,«tot, and Vhor = 



'v%. + Vy. We subtract the power maps of the quiet simulation computed in the same 
bandpasses to reduce the realization noise. The frequency filters used to recover the power maps and the 
azimuthally averaged power profiles (about the flux tube center) obtained subsequently are shown in Figure[5] 
Noteworthy aspects are that Whor contains the most intense halos, Vtot shows a dramatic increase in RMS 
power in the range 4-5 mHz around the 'umbral' region of the flux tube (defined as within a distance of 8 
Mm from the center of the flux tube), while v z displays limited shifts in the RMS in comparison to the rest. 
These effects (or some fraction thereof) could be attributed to changes in the eigenfunctions caused by the 
magnetic fields. The appropriate identification of the nature of these increments and decrements is evidently 
an important issue. 

Another set of power maps is displayed in Figure [H The upper set of panels contains the power maps of 
the 600 G flux tube whereas the lower two rows show the results from a simulation with a more realistically 
endowed sunspot: a 3000 G flux tu be (simulation size: 200 x 200 x 35 Mm 3 ). The flux tube configuration 
is very similar to that discussed in ICameron. Gizon. fc Daiffallahl (|2007l ) ; consequen tly, we do not s how i t 
here. The middle images are strikingly similar in structure to the observations of iMoretti et alJ (|2007l ). 
who see power increasing progressively with frequency (Figure 1 of their paper). We see a large decrease 
in the RMS power as felt by the pressure fluctuations (interpreted crudely as intensity) in the lowest set 
of panels. There is so me qualitative agreement between the simulations and the intensity observations 



More tti et al. ( 2007); h owever, the high resolution Hinode measurements of intensity in active regions 



Nagashima et al.l ( 2007) are unfortunate ly not so easily woven into this computational web. Intensity 



by 
by 

observations, as note lNagashima et alJ (|2007h . are far more difficult to interpret than those in velocity because 
of its sensitivity to the ionization, pressure, density etc., and the lack of a one to one correspondence with a 
simple thermodynamic variable. 

Acoustic halos around the edges of active regions have been wi dely observed (e .g. Braun et al. 19921 : 
Brown et all Il992l : balthasar et al.lll998t bonea et all I2OO0I ) . While iBalthasar et all |l998h have reported 
enhancements in oscillation velocity power also within magnetic regions and in low frequency (~ 2 mHz) 
bandpasses, a large number of other observations seem to show halos onl y in a high frequency ban dpass 
and in predominantly weakly magnetic areas surrounding the active region (jHindman fc Brownlll9 98). It is 
interesting to note that some qualitative features also seen in observations are reproduced in the simulations: 

(1) at the edge of the flux tube (at ~ 19 — 20 Mm in Figure [5j |B| ~ 7 — 12 G), only the highest frequency 
bandpass shows a faint power enhancement, of the order of 2 - 3 % in v z and even less in the other components, 

(2) the increase in wave power in the umbra of the flux tube in the 4-5 mHz bandpass (panel d of Figure [5]) 
is similar to enhancements seen in the magnetic cores of active regions ( Balthasar et al. 19981 ). and (3) the 
enhancements grow with frequency , as seen in the simulation of the 3000 G flux tube of Figure [6] and in 
observations bv lMoretti et ajj (|2007h . It would be rather ludicrous to make quantitative comparisons between 



- 11 - 



observations and the simulations because of the simplified nature of these calculations: the lack of radiative 
heat transfer, realistic wave mode damping, a penumbra, convection, unmodeled atmospheric magnetic fields 
etc. 

The speculation that enhanced seismic emission in t he vicinity of active regions may be the causative 
mecha nism of the acoustic halo goes back to the work of Brown et al. ( 19921 ). More recently, Donea et al 



(|2000[ ) have drawn similar conclusions from holography related analyses of active region observations. How- 
ever, this theory does not explain the wave power increase in the simulations because in our calculations, wave 
source amplitudes are statistically homogeneously distributed in space (in the horizontal directions) with 
the exception of areas close to the boundaries. Magnetic regions reconfigure the energy of the background 
medium. Therefore, the presence of sources in the interior of the flux tube essentially complicates matters 
because the incipient waves may have energies unlike waves in quiet regions. Moreover, the relative locations 
of the j3 = 1 layer with respect to the acoustic reflection zone, the r = 1 line, and the sources probably play 
an extremely important role in determining the wave energy distribution as a function of frequency. All the 
variables in the simulation are extracted at constant geometrical height (200 km above the photosphere), 
clearly a simplification incongruent with reality. Whether the observation height is a significant contributor 
is yet to be determined. Further investigations are currently in progress and will be the focus of a future 
paper. 



4. DISCUSSION 

We have discussed and validated a numerical method to systematically study linear MHD interactions 
in the context of helioseismology. The importance of including the ambient atmospheric magnetic field in 
the vicinity of magnetic flux concentrations is underlined here. Through a phenomenological model of the 
gradient smoothing that the ambient magnetic field presumably effects, we have shown that there can be 
significant differences in estimates of the oscillation velocity inside active regions. Thus, forward models 
that attempt to recover the magnetic field distribution based on shifts in travel times or other helioseismic 
metrics must in fact address this issue. Computational studies pertaining to oscillation power reduction in 
active regions are also quite sensitive to these effects. 

Results from simulations of waves interacting with 600 and 3000 G strong flux tubes are discussed in 
some detail. Not only is a significant reduction in wave power observed but a halo that surrounds the flux tube 
is also seen. Many features in the velocity observations of active regions are reproduce d by the simu l ations . 



High frequency wa ve p ower halos ar e also observed to envelope solar active regions; iBraun et al.l (|l992l ). 



Brown et al. (1992) and Donea et al. ( 2000h suggest enhanced seismic emission in the vicinity as being the 



causative mechanism. However, the simulations contain no such seismic enhancements, indicating that the 
physics behind the formation of the halo is possibly governed by MHD phenomena. A theory to explain the 
appearance of these excess oscillations will be discussed in a future publication. 

Using the techniques described here, we wish to develop helioseismically consistent forward models of 
thin flux tubes and sunspots. In the context of thin flux tube models, preli minary investigations have al - 
ready shown that the peak flux tube magnetic field strengths of about 80 G |Duvall. Birch, fc Gizonl 2006 ) 



as observed by the MDI instrument are too small by two orders of magnitude to cause the observed wave 
phase shifts. This is a consequence of the relatively low resolution of MDI, which is unable to capture the 
100-200 km sized flux tubes (Tom Duvall, Jr. 2007; Tom Bogdan 2007; Robert Cameron 2007, various private 
communications). Simulations with such small features can be computationally challenging due to rcsolu- 
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Fig. 5. — Azimuthal averages of the frequency filtered wave power maps (600G case). Panel a shows four 
filters with bandpasses 2 - 3, 3 - 4, 4 - 5, and 5-6 mHz. Panels b, c, and d display the azimuthally averaged 
(about the center of the flux tube) normalized noise-subtracted power maps of the quantities v z ,Vtot, and 
Each power difference is normalized by the mean value of the RMS of the correspondingly 
filtered quiet simulation. v z exhibits the least change in the RMS of all the variables shown here. Note that 
Whor and vtot are more difficult to interpret because sign information is lost (u?' as , ^to™ > 0). The two 
lines parallel to the y-axis in panels b,c, and d show the magnetic field strength and inclination at these 
locations. 



tion restrictions and the associated computational overhead. However, interesting sub-w avelength physics 
associated with thin flux tubes, namely the near-field evanescent modes (the jacket, e.g. iBogdan fc Callv 
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Fig. 6. — Changes in acoustic power caused due to MHD interactions as observed in v z , the vertical velocity 
(upper two rows), and p, the pressure. The color scale is fixed to the range [—50, +15]%, where each map 
has been normalized by the value of the quiet power in that frequency range. In both simulations, the 
source distributions are spatially uniform. The qualitative and quantitative differences seen between the 
two cases (images in the upper two rows) could be ascribed to the magnetic field strength, the location of 
(3=1 layer, and the source depth. We use pressure fluctuations as a proxy for observations in intensity. 
Although more careful studies are required to deconstruct these results into participat ory elements, it can 
be said that the compu tations reproduce many features in the velocity observations (e.g. lMoretti et al.ll2007 ; 



Nagashima et al. 2007). 



1995c iHanasoge et al.ll2007a) can be studied in greater detail with these simulations. These investigations 
are exciting, especially seen in the context of the availability of high quality observations and the upcoming 
Solar Dynamics Observatory (SDO) mission. 
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A. Validation: 2D analytical solution 



Take a 2D slab of finite thickness (L, L). Let the coordinates be labeled (x, z) and assume the presence a 
background magnetic field of the form Bo = Bo(x)e z . The background density is assumed to be unchanged 
by the magnetic field and is spatially non-varying; the pressure po is adjusted so that a pressure balance is 
achieved. We choose a velocity of the form, v = (v x e x + v z e z ) exp [i(kz — tot)], where v x = v x (x), v z = v z (x), 
k the wavenumber, li the frequency and t, time. Background quantities are denoted by the subscript 0. The 
magnetic field and pressure fluctuations are denoted by B and p respectively. Since this solution is used to 
validate the code, we use the linearized ideal MHD equations, which are equations ((T|) through (JSJ) without 
the boundary dissipative T or Lorentz force controller £ terms; we also set the source term S = 0. Starting 
with the adiabatic energy equation (upon incorporating the continuity equation), we have: 



dtp 
c 2 

dxPo 



Po 



-d x p 



-c 2 p Q V-v - v x d x p , 

Po 
-9 X 
1 



Po 



-d z p 



iujp 
d x v x d x 
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LOpo 



2 



Bl 



(d x V x + ikv z ) - Yip Q {dlv x + ikd x v z ) 

i[kz—ujt\ 



-TiPo(d x v x + ikv z ) + v x d x 



^i[kz— ujt\ 
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(A3) 



(A4) 
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where equation (|A3|) is the pressure distribution created by balancing the Lorentz force due to the background 
magnetic field and T\ is the first adiabatic index. Moving on to the x-momentum equation, and applying 
equation (|A4[) . 
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Similarly, upon the application of equation (I A5[) in the z-momentum equation, it may be verified that 



ikv z — 



2 2 k 2 



d x v x , 



leading to the relation 
where, 



d x v x + ikv z = r)(x)d x v x , 



r](x) 



1 



(A7) 
(A8) 
(A9) 



Upon further manipulation, a second-order differential equation for the eigenfunction v x may be obtained: 

d 2 x v x + Q(x)d x v x + <5>{x)v x = 0, (AlO) 



- 15 - 



where, 



e(x) 



Po^ 2 - 



k 2 Bl 



Bl 



(All) 
(A12) 



Tift)?? - 

Equation (|A10[) was solved using the MATLAB boundary value problem solver bvp4c. The boundary 
conditions were chosen to be u x U=o,£ = 0, with the additional condition d x v x \x=o = 1 required to solve for 
the eigenvalue u. Because of the linearity of the problem, there is no loss of generality due to this third 
condition. We show a sample eigenfunction calculation in Figure[7]for the resonant mode with v = 5.09 mHz; 
theory and simulation show good agreement. The background magnetic field was chosen to be -Bo = by/2x, 
b = 71.5 G Mm~ 1/2 , T x = 1.5, p = 1.21 x 10 5 -.Bg/2 dyne cm" 2 ,^ = 2.78 x 1(T 7 g cm" 3 , with x expressed 
in Mm. 
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Fig. 7. — The analytically computed (solid line) and numerically simulated (dot-dash line) eigenfunctions 
for v = 5.09 mHz, k — 0. For convenience, the background magnetic was chosen as Bq = b\/2x. Although 
not shown here, we have also tested the simulation at non-zero values of k and found good agreement. 
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